#核心论文代码验证
using Plots, ComplexPhasePortrait, ApproxFun, SpecialFunctions
gr();

u is in complex plane.

λ_{0}^{2} = (1 + u M_{0})^{2} - u^{2}

$\lambda_0^2 =(1 + uM_0)^2-u^2$

λ_{1}^{2} = C_{1}^{2} (1 + u M_{0})^{2} - u^{2}

$\lambda_1^2 =C_1^2(1 + uM_0)^2-u^2$ The principal branches of the square roots are used so that the branch cuts of

λ_{0}^{\pm}

$\lambda_0^\pm$ go from

u_{0}^{\pm}

$u_0^\pm$ to

\infty

$\infty$ and the branch cuts of

λ_{1}^{\pm}

$\lambda_1^\pm$ go from

u_{1}^{\pm}

$u_1^\pm$ to

\infty

$\infty$ where we have defined

u_{0}^{\pm} = 1 / (M_{0} \pm 1)

$u_0^\pm = 1/(M_0\pm1)$ and

u_{1}^{\pm} = C_{1} / (M_{1} C_{1} \pm 1)

$u_1^\pm = C_1/(M_1C_1\pm1)$ .

#u is in complex plane,for lanta_0,it has branch cuts
M0=0.3
phaseplot(-3..3, -3..3, z -> sqrt((1-z*M0)^2-z^2))

The fourier transform of the hard-wall boundary condtion on the centre body is simply $\partial \beta / \partial r\left ( h,u \right )=0$ 满足. To satisfy this boundary condition, we choose to write the solution of Bessel eqution as follows:

β (r, u) = {\begin{cases} A (u)) H_{m}^{(1)} (λ_{0} w r) & if r > 1 \\ B (u) [Y_{m}^{^{'}} (λ_{1} w h) J_{m} (λ_{1} w r) - J_{m}^{^{'}} (λ_{1} w h) Y_{m} (λ_{1} w r)] & if h < r < 1 \end{cases}

$\beta( r,u )=\begin{cases} A(u))H_m^{(1)}(\lambda_0wr) & \text{ if } r>1 \\ B(u)[Y_m^{'}(\lambda_1wh)J_m(\lambda_1wr)-J_m^{'}(\lambda_1wh)Y_m(\lambda_1wr)] & \text{ if } h<r<1 \end{cases}$

where $J_m$ and $Y_m$ are the Bessel and Neumann functions of order m, and $H_m^{(1)}=J_m+ i Y_m$ is the Hankel function of the first type.

参考：关于bessel的导数：https://www.boost.org/doc/libs/1_57_0/libs/math/doc/html/math_toolkit/bessel/bessel_derivatives.html http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/20/ShowAll.html

J_{m}^{^{'}} = \frac{J_{m - 1} - J_{m + 1}}{2} Y_{m}^{^{'}} = \frac{Y_{m - 1} - Y_{m + 1}}{2}

$J_m^{'}={J_{m-1}-J_{m+1}\over 2} \\ Y_m^{'}={Y_{m-1}-Y_{m+1}\over 2}$

#首先让我们大致看一下bessel函数的曲线形式
#与https://dlmf.nist.gov/10.3一致

𝐽𝑚=(m,x) -> besselj(m, x)
𝑌𝑚=(m,x) -> bessely(m, x)
x = range(0.01,stop=10,length=1000)
plot(x, 𝐽𝑚.(0,x);label="J_0",ylims = (-1,1))
plot!(x, 𝐽𝑚.(1,x);label="J_1")
plot!(x, 𝑌𝑚.(0,x);label="Y_0")
plot!(x, 𝑌𝑚.(1,x);label="Y_1")

m=1
𝐽′𝑚=(m,x) -> 1/2*(besselj(m-1, x)-besselj(m+1, x))
𝑌′𝑚=(m,x) -> 1/2*(bessely(m-1, x)-bessely(m+1, x))
x = range(0.01,stop=10,length=1000)
plot(x, 𝐽′𝑚.(m,x);label="J'_1")
plot!(x, 𝑌′𝑚.(m,x);label="Y'_1",ylims = (-1,1))

f = z -> hankelh1(3,z)
# set up plotting grid
z = range(0.01 ,stop=100,  length=1000)
plot(z, real(f.(z)))

#验证汉克函数
𝐻𝑚=(m,z) -> besselj(m, z)+im*bessely(m, z)
@show 𝐻𝑚(2, 1+im)
@show hankelh1(2,1+im)
phaseplot(-3..3, -3..3, z -> hankelh1(2,z))
#验证汉克函数在r趋近与无穷大时，叙述为正和负的区别
@show hankelh1(1,-100im)
@show hankelh1(1,+100im)

𝐻𝑚(2, 1 + im) = -0.5357570706365329 - 0.22597037902118686im
hankelh1(2, 1 + im) = -0.5357570706365328 - 0.22597037902118686im
hankelh1(1, -100im) = 1.30837515418464e26 - 2.1367387806763253e42im
hankelh1(1, 100im) = -2.9792874198947446e-45 - 1.8242874012570376e-61im

-2.9792874198947446e-45 - 1.8242874012570376e-61im

关于汉克函数的导数 https://keisan.casio.com/menu/system/000000000760

必须要指出的是，hankel函数用 $Jm+im*Ym$ 近似的方法只能针对与实数轴，对于过大的正虚数，imag》20im，则将不成立! 原因是，hankelh1(1, 20im)-->0,而besselj(1, 20im)-->∞&bessely(1, 20im)-->∞，因此造成很大的数值误差. 举例说明：

#验证hankel函数的导数
𝐻1𝑚=(m,z) -> hankelh1(m, z)
H1𝑚=(m,z) -> besselj(m, z)+im*bessely(m, z)
𝐻′1𝑚1=(m,z) ->hankelh1(m-1, z)-m/z*hankelh1(m, z)
H′1𝑚1=(m,z) ->besselj(m-1, z)+im*bessely(m-1, z)-m/z*(besselj(m, z)+im*bessely(m, z))
#H′1𝑚2=(m,z) ->-(besselj(m+1, z)+im*bessely(m+1, z))+m/z*(besselj(m, z)+im*bessely(m, z))
#𝐻′1𝑚2=(m,z) ->-hankelh1(m+1, z)+m/z*hankelh1(m, z)

@show besselj(1,2+20*im)
@show bessely(1,2+20*im)
@show 𝐻1𝑚(1,2+20*im)
@show H1𝑚(1,2+20*im)
@show besselj(1,2-20*im)
@show bessely(1,2-20*im)
@show 𝐻1𝑚(1,2-20*im)
@show H1𝑚(1,2-20*im)

zr=range(-40,stop=40,length=1000)
zi=ones(1000,1)*im*20
z=zr+zi
r_j = z -> real(besselj(1,z))
i_j = z -> imag(besselj(1,z))
r_y = z -> real(bessely(1,z))
i_y = z -> imag(bessely(1,z))
r_h=z -> real(hankelh1(1,z))
i_h=z -> imag(hankelh1(1,z))
#plot(zr,r_j.(z);label="r-j")
#plot!(zr,i_j.(z);label="i-j")
#plot!(zr,r_y.(z);label="r-y")
#plot!(zr,i_y.(z);label="i-y")
plot(zr,r_h.(z);label="T-r-h")
plot!(zr,i_h.(z);label="T-i-h")

#plot!(zr,r_j.(z)+i_y.(z);label="F-r-h")
#plot!(zr,i_j.(z)+r_y.(z);label="F-i-h")

besselj(1, 2 + 20im) = 3.931553698919195e7 - 1.576352471413e7im
bessely(1, 2 + 20im) = 1.5763524714130001e7 + 3.931553698919195e7im
𝐻1𝑚(1, 2 + 20im) = 1.727569618545747e-10 - 3.311815065025103e-10im
H1𝑚(1, 2 + 20im) = 0.0 + 1.862645149230957e-9im
besselj(1, 2 - 20im) = 3.931553698919195e7 + 1.576352471413e7im
bessely(1, 2 - 20im) = 1.5763524714130001e7 - 3.931553698919195e7im
𝐻1𝑚(1, 2 - 20im) = 7.86310739783839e7 + 3.1527049428260002e7im
H1𝑚(1, 2 - 20im) = 7.86310739783839e7 + 3.152704942826e7im

M0=0.3
m=2
w=1+im
λ0=(u)->sqrt((1-u*M0)^2-u^2)
𝐻1𝑚=(u) -> besselj(m,u)+im*bessely(m,u)
𝐻2𝑚=(u) -> hankelh1(m,u)
𝐻1𝑚_λ0w=(u) -> besselj(m, λ0(u)*w)+im*bessely(m, λ0(u)*w)
𝐻2𝑚_λ0w=(u) -> hankelh1(m, λ0(u)*w)
𝐻′1𝑚1_λ0w=(u) ->hankelh1(m-1, λ0(u)*w)-m/(λ0(u)*w)*hankelh1(m, λ0(u)*w)
H′1𝑚1_λ0w=(u) ->besselj(m-1, λ0(u)*w)+im*bessely(m-1, λ0(u)*w)-m/(λ0(u)*w)*(besselj(m, λ0(u)*w)+im*bessely(m, λ0(u)*w))

phaseplot(-20..20, -40..40, 𝐻′1𝑚1_λ0w)
phaseplot(-20..20, -40..40, H′1𝑚1_λ0w)

m=2
w=1+im
λ0=(u)->sqrt((1-u*M0)^2-u^2)
𝐻1𝑚_λ0w=(u) -> besselj(m, λ0(u)*w)+im*bessely(m, λ0(u)*w)
𝐻2𝑚_λ0w=(u) -> hankelh1(m, λ0(u)*w)
phaseplot(-20..20, -20..20, 𝐻2𝑚_λ0w)
phaseplot!(-20..20, -20..20, 𝐻1𝑚_λ0w)

The fourier transform of the hard-wall boundary condtion on the centre body is simply $\partial \beta / \partial r\left ( h,u \right )=0$ 满足. To satisfy this boundary condition, we choose to write the solution of Bessel eqution as follows:

β (r, u) = {\begin{cases} A (u)) H_{m}^{(1)} (λ_{0} w r) & if r > 1 \\ B (u) [Y_{m}^{^{'}} (λ_{1} w h) J_{m} (λ_{1} w r) - J_{m}^{^{'}} (λ_{1} w h) Y_{m} (λ_{1} w r)] & if h < r < 1 \end{cases}

$\beta( r,u )=\begin{cases} A(u))H_m^{(1)}(\lambda_0wr) & \text{ if } r>1 \\ B(u)[Y_m^{'}(\lambda_1wh)J_m(\lambda_1wr)-J_m^{'}(\lambda_1wh)Y_m(\lambda_1wr)] & \text{ if } h<r<1 \end{cases}$

where $J_m$ and $Y_m$ are the Bessel and Neumann functions of order m, and $H_m^{(1)}=J_m+ i Y_m$ is the Hankel function of the first type.

那么问题转化为，在复平面范围内， $\lambda$ 避开branch cut的情况下，如何保证 $\lambda_0wr$ 虚数部分始终为正数？其中， $w=w_r+iw_i=|w|e^{i\varepsilon }, 0\leqslant \varepsilon \leqslant \pi/2$ ,即 $w_r,w_i$ 均大于0

θ=range(0.0,stop=π/2-0.01,length=100)
plot(θ,tan.(θ))

wr = 0:1:100
wi = 0:1:100
f = z -> z
r = z -> abs(f(z))
θ = z -> angle(f(z))
plot( contourf(wr, wi, r.(wr' .+ im.*wi); title="abs"),
      contourf(wr, wi, θ.(wr' .+ im.*wi); title="phase"))

M0=0 #M0=0,u±=±1
wr = 4
wi = 0
w=wr+ im*wi
phaseplot(-3..3, -3..3, z ->w*sqrt((1-z*M0)^2-z^2))

#u is in complex plane,for lanta_0,it has branch cuts
#λ0=sqrt((1-z*M0)^2-z^2)
#λ1=sqrt((C1-z*M0)^2-z^2)
M0=0 #M0=0,u±=±1
M1=0.5
C1=0.8
wr = 4
wi = 0
w=wr+ im*wi

xx = -10:0.01:10
yy = -10:0.01:10
f = z -> w*sqrt((1-z*M0)^2-z^2)

r = z -> abs(f(z))
θ = z -> angle(f(z))
image=z -> imag(f(z))
plot( contourf(xx, yy, r.(xx' .+ im.*yy); title="abs"),
      contourf(xx, yy, θ.(xx' .+ im.*yy); title="phase"),
      contourf(xx, yy, image.(xx' .+ im.*yy); title="image"))

#u is in complex plane,for lanta_0,it has branch cuts
#λ0=sqrt((1-z*M0)^2-z^2)
#λ1=sqrt((C1-z*M0)^2-z^2)
M0=0.3 #M0=0,u±=±1
M1=0.5
C1=0.8
wr = 4
wi = 4
w=wr+ im*wi

xx = -10:0.01:10
yy = -10:0.01:10
f = z -> w*sqrt((1-z*M0)^2-z^2)
#r = z -> abs(f(z))
#θ = z -> angle(f(z))
image=z -> imag(f(z)) #求f函数的虚数部分

tanε=wi/wr
u1ip=range(10,stop=-10,length=10)
u1rp=u1ip./-tanε+ones(10,1)*1/(M0+1)
u1im=range(10,stop=-10,length=10)
u1rm=-u1im./tanε+ones(10,1)*1/(M0-1)
u2ip=range(10,stop=-10,length=10)
u2rp=u2ip./-tanε+ones(10,1)*C1/(C1*M0+1)
u2im=range(10,stop=-10,length=10)
u2rm=-u2im./tanε+ones(10,1)*1/(C1*M0-1)
plot(contourf(xx, yy, image.(xx' .+ im.*yy); title="image"))
plot!(u1rp,u1ip,aspect_ratio=:equal;label="R1+0",arrow=1)
plot!(u1rm,u1im,aspect_ratio=:equal;label="R1-0",arrow=1)

#plot!(u2rp,u2ip,aspect_ratio=:equal;label="R2+0")
#plot!(u2rm,u2im,aspect_ratio=:equal;label="R2-0")


#phaseplot(-3..3, -3..3, z ->w*sqrt((1-z*M0)^2-z^2))

Finally,we can derive the Winer-Hopf equation:

G (u) = G_{-} (u) + G_{+} (u) = w \frac{K (u) F_{+} (u)}{λ_{0} λ_{1}}

$G(u)=G_-(u)+G_+(u)=w{K(u)F_+(u) \over \lambda_0 \lambda_1}$ where K is the Wiener-Hopf kernel:

K (u) = D_{1} (1 - u M_{1})^{2} λ_{0} \frac{Y_{m}^{^{'}} (λ_{1} w h) J_{m} (λ_{1} w) - J_{m}^{^{'}} (λ_{1} w h) Y_{m} (λ_{1} w)}{Y_{m}^{^{'}} (λ_{1} w h) J_{m}^{^{'}} (λ_{1} w) - J_{m}^{^{'}} (λ_{1} w h) Y_{m}^{^{'}} (λ_{1} w)} - (1 - u M_{0})^{2} λ_{1} \frac{H_{m}^{^{'}} (λ_{0} w)}{H_{m}^{(1)^{'}} (λ_{0} w)}

$K(u)=D_1(1-uM_1)^2 \lambda_0 {Y_m^{'}(\lambda_1wh)J_m(\lambda_1w)-J_m^{'}(\lambda_1wh)Y_m(\lambda_1w)\over Y_m^{'}(\lambda_1wh)J_m^{'}(\lambda_1w)-J_m^{'}(\lambda_1wh)Y_m^{'}(\lambda_1w)} -(1-uM_0)^2\lambda_1{H_m^{'}(\lambda_0w)\over H_m^{(1)'}(\lambda_0w)}$ 接下来，对

K (u)

$K(u)$ 的zeros和poles进行分析。

M0=0.25
M1=0.5
D1=C1=1
m=2
w=4+4im
h=0.1
λ0=(u)->sqrt((1-u*M0)^2-u^2)
λ1=(u)->sqrt((C1-u*M1)^2-u^2)
tanε=imag(w)/real(w)
u1ip=range(10,stop=-10,length=10)
u1rp=u1ip./-tanε+ones(10,1)*1/(M0+1)
u1im=range(10,stop=-10,length=10)
u1rm=-u1im./tanε+ones(10,1)*1/(M0-1)
u2ip=range(10,stop=-10,length=10)
u2rp=u2ip./-tanε+ones(10,1)*C1/(C1*M1+1)
u2im=range(10,stop=-10,length=10)
u2rm=-u2im./tanε+ones(10,1)*1/(C1*M1-1)

#####
𝐽𝑚_λ1w=(u) -> besselj(m, λ1(u)*w)
𝑌𝑚_λ1w=(u) -> bessely(m, λ1(u)*w)
#𝐽𝑚_λ0w=(u) -> besselj(m, λ0(u)*w)
#𝑌𝑚_λ0w=(u) -> bessely(m, λ0(u)*w)
𝐽′𝑚_λ1w=(u) -> 1/2*(besselj(m-1, λ1(u)*w)-besselj(m+1, λ1(u)*w))
𝑌′𝑚_λ1w=(u) -> 1/2*(bessely(m-1, λ1(u)*w)-bessely(m+1, λ1(u)*w))
𝐽′𝑚_λ1wh=(u) -> 1/2*(besselj(m-1, λ1(u)*w*h)-besselj(m+1, λ1(u)*w*h))
𝑌′𝑚_λ1wh=(u) -> 1/2*(bessely(m-1, λ1(u)*w*h)-bessely(m+1, λ1(u)*w*h))
𝐻1𝑚1_λ0w=(u) -> hankelh1(m, λ0(u)*w)
𝐻′1𝑚1_λ0w=(u) ->hankelh1(m-1, λ0(u)*w)-m/(λ0(u)*w)*hankelh1(m, λ0(u)*w)
#𝐻′1𝑚2_λ0w=(u) ->-hankelh1(m+1, λ0(u)*w)+m/(λ0(u)*w)*hankelh1(m, λ0(u)*w)
K=(u)->D1*(1-u*M1)^2*λ0(u)*(𝑌′𝑚_λ1wh(u)*𝐽𝑚_λ1w(u)-𝐽′𝑚_λ1wh(u)*𝑌𝑚_λ1w(u))/(𝑌′𝑚_λ1wh(u)*𝐽′𝑚_λ1w(u)-𝐽′𝑚_λ1wh(u)*𝑌′𝑚_λ1w(u))-(1-u*M0)^2*λ1(u)*𝐻1𝑚1_λ0w(u)/𝐻′1𝑚1_λ0w(u)
phaseplot(-15..15,-10..10, K)
plot!(u1rp,u1ip,aspect_ratio=:equal;label="R1+0",arrow=1,legend=:bottomright)
plot!(u1rm,u1im,aspect_ratio=:equal;label="R1-0",arrow=1)
plot!(u2rp,u2ip,aspect_ratio=:equal;label="R2+0",arrow=1)
plot!(u2rm,u2im,aspect_ratio=:equal;label="R2-0",arrow=1)